A New Aryl Hydrocarbon Receptor Homology Model Targeted to Improve Docking Reliability

The aryl hydrocarbon receptor (AhR) is a ligand-dependent, basic helix-loop-helix Per-ARNTSim (PAS) containing transcription factor that can bind and be activated by structurally diverse chemicals, including the toxic environmental contaminant 2,3,7,8-tetrachlorodibenzo-p-dioxin (TCDD). As no experimentally determined structures of the AhR ligand binding domain (LBD) are available and previous homology models were only derived from apo template structures, we developed a new model based on holo X-ray structures of the hypoxia-inducible factor 2α (HIF-2α) PAS B domain, targeted to improve the accuracy of the binding site for molecular docking applications. We experimentally confirmed the ability of two HIF-2α crystallographic ligands to bind to the mAhR with relatively high affinity and demonstrated that they are AhR agonists, thus justifying the use of the holo HIF-2α structures as templates. A specific modeling/ docking approach was proposed to predict the binding modes of AhR ligands in the modeled LBD. It was validated by comparison of the calculated and the experimental binding affinities of active THS ligands and TCDD for the mAhR and by functional activity analysis using several mAhR mutants generated on the basis of the modeling results. Finally the ability of the proposed approach to reproduce the different affinities of TCDD for AhRs of different species was confirmed, and a first test of its reliability in virtual screening is carried out by analyzing the correlation between the calculated and experimental binding affinities of a set of 14 PCDDs. Introduction The aryl hydrocarbon receptor (AhR)a is a basic helix-loop-helix (bHLH), Per-ARNT-Sim (PAS) containing ligand-dependent transcription factor that induces the expression of a large battery of genes and produces diverse biological and toxic effects in a wide range of species and tissues.1–4 The best-characterized high affinity ligands include a variety of toxic halogenated aromatic hydrocarbons (HAHs), such as the polychlo-dibenzo-p-dioxins (PCDDs), dibenzofurans (PCDFs), and biphenyls (PCBs), and numerous polycyclic *To whom correspondence should be addressed. Phone: (+39)0264482821. Fax: (+39)0264482839. [email protected]. Supporting Information. Four additional Figures are provided as Supporting Information (SI): Figure SI1. Alignment of the LBD domains of mouse, human and rat AhRs against HIF-2α. Figure SI2. Cα-trace representation of the 100 mAhR homology models generated by MODELLER. Figure SI3. Results of ensemble docking of THS-017 and THS-020 in the mAhR homology models. Figure SI4. Results of PCDD ensemble docking in the rtAhR LBD homology models. This information is available free of charge via the Internet at http://pubs.acs.org/. NIH Public Access Author Manuscript J Chem Inf Model. Author manuscript; available in PMC 2012 November 28. Published in final edited form as: J Chem Inf Model. 2011 November 28; 51(11): 2868–2881. doi:10.1021/ci2001617. N IH PA Athor M anscript N IH PA Athor M anscript N IH PA Athor M anscript aromatic hydrocarbons (PAHs) and PAH-like chemicals,5–7 all widespread classes of environmental contaminants. Moreover, a number of natural, endogenous and synthetic AhR agonists and antagonists whose structure and physicochemical characteristics are different from those of the prototypical HAH and PAH ligands have been identified as lower affinity ligands and moderately potent inducers of AhR-dependent gene expression.7–11 Among the various protein domains responsible for the AhR functional activities, PAS B (one of the two structural repeats in the PAS domain) is the one responsible for ligand binding and it is also involved in binding to the chaperone heat shock protein 90 (hsp90).12 Although the AhR signal transduction pathway has been studied for many years,1–4 several unanswered questions remain. Major issues are the actual spectrum of ligands, how they can bind to the AhR, and how ligand binding to the ligand binding domain (LBD) results in activation of the AhR and AhR-dependent gene expression. A molecular understanding of these events would require detailed structural information about the AhR PAS B LBD. However, neither X-ray nor NMR structures of the bound or unbound AhR have been determined to date. Since the first crystal structures of distant homologous proteins belonging to the PAS superfamily became available, we started developing theoretical models of the AhR LBD by homology modeling techniques and the results provided an initial framework to make hypotheses on LBD characteristics and the mechanisms of AhR functionality.9,13,14 The latest model of the mouse AhR (mAhR) LBD we proposed14 was built using the NMR structures of the PAS B domains of the human hypoxia-inducible factor 2α (HIF-2α)15 and of the human ARNT,16 both in the apo form, as templates. That was the most reliable among the models developed in our group, since the template domains show the highest degree of sequence identity and similarity with the AhR PAS B among all the PAS structures reported to date. Moreover, the full length template proteins are members of the bHLH/PAS family of transcriptional factors and functionally related to the AhR.4 On the basis of model-driven site-directed mutagenesis and AhR functional analysis, the buried cavity in the core of the domain was confirmed as the site involved in ligand binding.14 Moreover, analysis of the LBD models of several mammalian AhRs indicated that the physico-chemical characteristics of the binding cavities are remarkably conserved in all AhRs with high affinity for the AhR ligand 2,3,7,8-tetrachlorodibenzo-p-dioxin (TCDD).17 Mutagenesis of the conserved residues, followed by AhR functional analysis, allowed identification of the “TCDD-binding fingerprint”, i.e. the group of residues necessary for optimal TCDD binding.17 On the basis of our findings, other authors recently reported the development of AhR LBD homology models based on the same apo structure of HIF-2α and proposed the use of these models for molecular docking applications.18–22 Although the use of homology models of the target protein, in addition to experimental structures, has greatly extended the applicability of molecular docking approaches,23–25 the use of good quality models, in particular of the binding site, is crucial for the prediction of reliable binding poses.24, 26–28 In fact, it was proven that the success of docking calculation decreases with the quality of the receptor structure and, in particular, with the decreasing ability of the structure to aAbbreviations: AhR, aryl hydrocarbon receptor; mAhR, Mus Musculus (mouse) aryl hydrocarbon receptor; huAhR, Homo Sapiens (human) aryl hydrocarbon receptor; rtAhR, Rattus Norvegicus (rat) aryl hydrocarbon receptor; ARNT, aryl hydrocarbon receptor nuclear translocator; HIF-2α, hypoxia-inducible factor 2α; hsp90, heat shock protein 90; bHLH, basic Helix-Loop-Helix; DRE, dioxin responsive element; LBD, ligand binding domain; PAS, Per-ARNT-Sim; TCDD, 2,3,7,8-tetrachlorodibenzo-p-dioxin; HAH, halogenated aromatic hydrocarbons; PCDD, polychlorodibenzo-p-dioxin (MCDD, monochloro; DCDD, dichloro; TrCDD, trichloro; TCDD, tetrachloro; PeCDD, pentachloro; Hx, exachloro; O, octachloro); PCDF, polychlorodibenzofuran; PCB, polychlorobiphenyl; PAH, polycyclic aromatic hydrocarbons; PDB, Protein Data Bank; NMR, Nuclear Magnetic Resonance; RMSD, root mean square deviation; HM, homology model. Motto et al. Page 2 J Chem Inf Model. Author manuscript; available in PMC 2012 November 28. N IH PA Athor M anscript N IH PA Athor M anscript N IH PA Athor M anscript reproduce the features of the active site that are important for ligand recognition. Modeled structures, even those having high sequence identity to their template structure, can have improperly placed side chains in the binding site. This can be partially overcome by the use of holo template structures, that take into account the induced fit effects associated with the presence of a ligand in the binding cavity, thus greatly improving the quality of the modeled cavity and the reliability of the predicted binding modes of similar ligands.26, 29 Since all of the AhR LBD models developed to date were derived from apo template structures, their use in molecular docking could have some limitations. The recent depositions of new X-ray structures of the human HIF-2α:ARNT PAS B heterodimer co-crystallized with artificial ligands30,31 provided us with an opportunity to improve the quality of the modeled AhR binding site. With their studies, Gardner and coworkers demonstrated that other PAS domains, in addition to the AhR, are able to bind small-molecule ligands.30–32 Moreover, they highlighted that the crystal structures of the ligand-bound HIF-2α PAS B domain show larger binding cavities and have important conformational differences compared to the apo structures previously determined by NMR.30,31 Demonstration of the ability of these HIF-2α ligands to bind to the AhR LBD would suggest similarities in features of the two binding cavities that are involved in ligand recognition. This experimental finding would also justify the use of the new ligand-bound HIF-2α threedimensional structures to develop a new homology model of the AhR LBD, specifically targeted to improve the accuracy of the binding site for molecular docking applications aimed to reliably analyze the binding modes of the wide spectrum of chemicals that bind to the AhR. Here we demonstrate the binding and agonist activity of HIF-2α crystallographic ligands30,31 to the AhR, confirming the choice of using the above mentioned ligand-bound HIF-2α complexes as templates for the AhR LBD modeling. Using this improved model, we propose to predict the binding modes of AhR ligands in the modeled LBD by a specifically developed docking approach. First, validation of the homology modeling/docking approach is performed by comparison of the calculated and the experimental binding affinities of active THS ligands and TCDD for the mAhR. Moreover, the predicted TCDD binding mode within the binding cavity is confirmed by functional activity analysis using several mAhR mutants generated on the basis of the modeling results. Ligand binding to new homology models of the human AhR (huAhR) and the rat AhR (rtAhR) is also analyzed. With this information, the ability of the proposed approach to reproduce the different affinities of TCDD for AhRs of different species33 was verified and a first test of its reliability in virtual screening applications is carried out by analyzing the correlation between the calculated and experimental binding affinities of a set of 14 PCDDs for the rtAhR.34,35 Results and Discussion THS ligands can bind to and activate the mAhR LBD A series of synthetic compounds that could bind to and be co-crystallized with a recombinant HIF-2α PAS B fragment31 has recently been reported. Considering that both HIF-2α and AhR proteins belong to the same superfamily and have PAS B domains with relatively high sequence identity and similarity, we hypothesized that these compounds might also bind to the AhR. Accordingly, we examined the ability of three of these compounds (THS-017, THS-020 and THS-044 (Figure 1)) to compete with [3H]TCDD for binding to in vitro synthesized mAhR.36 Preliminary experiments demonstrated that THS-017 and THS-020, but not THS-044, could competitively inhibit the specific binding of [3H]TCDD to the AhR (data not shown). Accordingly, to estimate the relative ligand Motto et al. Page 3 J Chem Inf Model. Author manuscript; available in PMC 2012 November 28. N IH PA Athor M anscript N IH PA Athor M anscript N IH PA Athor M anscript binding affinity of these new AhR ligands, competitive hydroxyapatite binding assays with in vitro synthesized mAhR was carried out with 2 nM [3H]TCDD in absence or presence of increasing concentrations of THS-017 and THS-020.36 These experiments revealed that THS-017 and THS-020 could competitively displace [3H]TCDD specific binding to the AhR in a concentration-dependent manner (Figure 2A). Relative binding affinity values (IC50) for THS-017 (~0.6 μM) and THS-020 (~4.2 μM) were determined by non-linear regression analysis of the competitive binding curves and demonstrate that these ligands have significantly lower affinity than that of TCDD (1 nM).37,38 While binding analysis confirms the ability of the compounds to directly interact with the AhR, it does not provide any information as to the functional activity of the ligand. To determine whether these ligands were AhR agonists, we examined their ability to stimulate DNA binding of in vitro synthesized mAhR (using gel retardation analysis) and AhR-dependent gene expression in cells in culture (using transient transfection). Both THS-017 and THS-020 could stimulate mAhR DNA binding in a concentration-dependent manner (Figure 2B) and the amount of induced THS-017/THS-020:AhR:ARNT:DNA complex formed at each concentration from multiple experiments was determined and expressed relative to that produced a maximal activating concentration of TCDD (Figure 2C). The relative binding affinity of THS-017 and THS-020 (Figure 2A) compares favourably with the relative potency of each to stimulate mAhR DNA binding (Figure 2B/C), with THS-017 having a slightly greater affinity and potency. The ability of these compounds to stimulate AhR transformation and DNA binding is consistent with their ability to act as AhR agonists. To confirm this, we examined the ability of THS-017 and THS-020 to stimulate AhR-dependent luciferase reporter gene expression. Cos-1 cells that had been transiently co-transfected with a mAhR expression vector and an AhR-responsive luciferase plasmid (pGudLuc6.1)39 were incubated for 24 hours with THS-017, THS-020 or several known AhR agonists (TCDD, 3methylcholanthrene (3MC), indirubin or alpha-naphthoflavone (aNF)) and luciferase activity determined. These results clearly demonstrate the ability of both THS-017 and THS-020 to stimulate AhR-dependent gene expression and although they were less potent than the well established AhR agonists TCDD, 3MC and indirubin, they were comparable to the lesser potent AhR agonist aNF. Together, these results demonstrate that THS-017 and THS-020 are AhR agonists with affinities and potencies comparable to many well-characterized AhR ligands40 but they provide strong justification for use of the crystal structures of THS-bound HIF-2α as improved templates for modeling the AhR LBD. mAhR homology modeling Homology modeling techniques were applied to predict the LBD structure (residues 278– 384) of the mAhR that takes into account the presence of ligands in its binding site. In fact, an improved description of the binding cavity in the bound form makes the mAhR model more suitable for accurate docking calculations and virtual screening applications26 than the models based on previously proposed apo template structures.14,17 This new homology model was derived by mono-template sequence alignment with the HIF-2α PAS B domain, for which multi-structural information exists. The availability of experimental structures of the HIF-2α PAS B domain co-crystallized with THS ligands (PDB ID: 3F1O,30 3H7W,31 3H8231) provides optimal template structures to improve the AhR PAS B homology model for docking applications. Their properties are summarized in Table 1. Despite our experiments demonstrated that only THS-017 and THS-020 act as AhR agonists, all the three template structures were used in the homology modeling. In fact, while the ligand binding site is highly conserved among the three crystallized THS:HIF-2α PAS B, the combined use of the three structures does provide a better statistical representation of the possible conformational variability of the whole protein in the bound form. Motto et al. Page 4 J Chem Inf Model. Author manuscript; available in PMC 2012 November 28. N IH PA Athor M anscript N IH PA Athor M anscript N IH PA Athor M anscript The sequence alignment of the mAhR and the HIF-2α template previously proposed14,17 was employed for modeling. Analysis of the alignment (Figure SI1) demonstrates that the major part of the AhR LBD can be modeled on the basis of all the three X-ray depositions of the template. A short region, corresponding to the Hβ/Iβ loop (the secondary structures attributed by DSSPcont41 are shown in the Figure), is solved only in one template structure (PDB ID: 3H82). The crystallographic structures of the three templates were prepared as described in the Experimental Section. One-hundred three-dimensional models for the mAhR LBD were generated using MODELLER42–45 and energetically ranked by the DOPE score46 (see Experimental Section). The template-bound structures of the THS ligands (Figure 1) were transferred into the mAhR homology models for the purpose of obtaining a binding cavity that takes into account ligand induced-fit effects. The first visual inspection of the resulting one-hundred mAhR homology models highlighted some conformational variability in the Hβ/Iβ loop and around the inserted Gly 313 (Figure SI2). The flexibility of the Hβ/Iβ loop can be ignored as it doesn’t influence the dimension of the binding site. In contrast, the conformational variability observed around Gly 313, that lines the ligand binding cavity, was further analyzed since the dynamic behavior of the binding site must be taken into account for docking applications. Based on these observations, a workflow for the selection of representative conformations of the mAhR binding site was established and the implemented selection funnel is presented in Scheme 1. Five conformational clusters for the AhR were identified and the almost regular membership distribution in clusters 1–4 confirmed a dynamic binding site behavior. Since only two conformers populate cluster 5, it can be considered statistically less representative in the description of the mAhR conformational variability and thus it was not included in subsequent modeling analyses and applications. PROCHECK47 validation of the selected representative conformations of each cluster indicated a good stereochemical quality, with 87 – 91% of residues belonging to the most favored areas of the Ramachandran plot and overall G-factors ranging from −0.26 to −0.19 (this index ranges from −0.5 to 0.3 for structures solved at 1.5 Å resolution). Moreover the ProSA z-scores48,49 values were between −4.14 and −3.87, within the range of values typical for native protein structures of similar size. Comparison of the apo and holo homology models Comparative analysis of the new mAhR homology models (holo mAhR) and the previously published model (apo mAhR14,17) was carried out to elucidate the effects of using holo template structures and including ligands during the binding site modeling. The global and the binding site RMSD of the holo mAhR representative conformations (in the following also named representative homology models, HM) were calculated relative to the apo mAhR model. The resulting values (Table 2) revealed that significant contributions to the global RMSD appear to derive from differences in the conformations of the binding site residues. Focusing on the comparative analysis of residues within 5 Å of the THS ligands in the holo model, several important differences between the apo and holo AhR LBD models appear. Particularly evident is the different conformation of the His 320 residue, whose side chain is projected into the binding site in the apo model while it is projected outside in the holo models. This directly results from the different structural depositions used in the homology modeling. In fact, in the sequence alignment (Figure SI1), AhR His 320 corresponds to His 282 in HIF-2α and it points away from the binding site in the utilized templates, while in the NMR apo HIF-2α structure (PDB ID: 1P9715), used in the previous homology modeling studies,14,17 His 282 points inward. However, it is conceivable that this conformational Motto et al. Page 5 J Chem Inf Model. Author manuscript; available in PMC 2012 November 28. N IH PA Athor M anscript N IH PA Athor M anscript N IH PA Athor M anscript effect is not related to a structural rearrangement of the binding site that occurs upon ligand binding but to crystallographic packing, since His 282 adopts the same outward conformation in all the available X-ray structures of HIF-2α, including the apo ones.30,31 Other conformational differences between the apo and holo template structures (i.e. side chain and backbone shifts) are also observed for residues lining the binding pocket, and they are clearly reflected in the respective mAhR homology models. In fact, the RMSD values between the residues identified as the “TCDD-binding fingerprint”17 in the holo and in the apo mAhR models (Table 2) indicate that significant conformational changes were induced in these side-chains by the use of holo template structures. Since the internal cavity volume is determined by the conformations assumed by the internal binding site residues, the CASTp50 web server51 was used to calculate and compare the accessible Connoly’s volumes52 for the main cavity of the HIF-2α templates, the published apo and the representative holo mAhR homology models. The results in Table 3 indicate a clear ligand induced-fit effect on the binding site extension, with inclusion of THS ligands in mAhR binding site almost doubling the cavity volume. Comparison of selected holo mAhR homology models with the HIF-2α templates reveals a more extended “channel-like” binding cavity. Superimposition of these models with the templates shows an outward shift of the Aβ strand that could explain the observed longer extension of the mAhR ligand binding site. Evaluation of the THS ligands binding mode in mAhR THS compounds were rigidly copied from the HIF-2α structural templates into the mAhR binding site during the homology modeling process. The representative conformations complexed with the HIF-2α ligands were subjected to energy minimization as described in the Experimental Section. This permitted relaxation of THS ligands in their new protein environment, in order to facilitate the evaluation of their possible binding mode in the mAhR LBD. All the energy minimized complexes were ranked using Glide extra precision (XP) scoring function53 and the representative complexes with the best scores (mAhR_HM79+THS-017, mAhR_HM70+THS-020, and mAhR_HM01+THS-044) were chosen for more focused structural analysis. An intramolecular hydrogen bond between the ligand aniline and nitro group, similar to the one described for the crystallographic HIF-2α/THS complexes,30,31 was observed also in the minimized complexes (Figure 3). This bond provides a high internal stabilization and favors a stable bent conformation of the ligands. Polar interactions of the amino and nitro groups with the His 285 and the Gln377 side chains anchor the central part of the ligand molecules. The trifluoromethyl group, that lies at the entrance of the cavity in a hydrophobic environment, is also subjected to some stabilizing polar interactions with Ser 330 and Ser 359, whereas the aryl ring is packed with the side chains of Phe 289 and Phe 345. As highlighted above, the mAhR LBD shows a more extended binding cavity than HIF-2α. Given the low calculated RMSD between the ligand structures preand post-energy minimization (Table 4), the binding conformations of the THS ligands appear almost conserved after the energy minimization. However, as both THS-017 and THS-020 contain a flexible spacer (methylene group), they can assume slightly more extended conformations in mAhR than in HIF-2α by projecting their heterocyclic rings deeper into the mAhR cavity, in a highly hydrophobic region. As shown in Figure 3, the chemical properties and the binding conformations observed for the heterocycles of THS-017, THS-020 and THS-044 in this zone seem to correlate with the relative order of their experimental affinities. In fact, the thiophene ring in THS-017 is highly hydrophobic and thus would be particularly well stabilized in this environment by favorable Van der Waals contacts. The furan group in THS-020 would also prefer hydrophobic surroundings but the harder oxygen atom introduces some system polarity; in the derived model only the Motto et al. Page 6 J Chem Inf Model. Author manuscript; available in PMC 2012 November 28. N IH PA Athor M anscript N IH PA Athor M anscript N IH PA Athor M anscript side chains of Thr 283, His 285 and Gln 377 could help in stabilizing this slightly polar system. In contrast, the placement of the polar and basic morpholine group of THS-044 would not be favored in this binding site area and this destabilization could explain its extremely low affinity for the AhR compared to the other THS ligands. Considering the experimental affinity data we obtained for the THS compounds in mAhR, THS-044, that does not bind mAhR, was excluded from further modeling investigations. Ensemble docking approach to study ligand binding to the AhR As discussed above, in mAhR homology models some conformational variability is observed around Gly 313. This residue and its surroundings line the ligand binding pocket and its dynamic behavior could influence the dimensions of the cavity and consequently ligand placement. For this reason, its dynamics should be taken into account during docking calculations. However, incorporating protein flexibility during molecular docking still remains one of the major challenges of this method.54,55 Although numerous programs allow induced-fit docking, these calculations still remain computationally expensive, time consuming and not suitable for virtual screening. Accordingly, ensemble docking56 was selected as the optimal approach to perform the “dynamic” docking of ligands in the mAhR homology model, using the representative homology models selected by cluster analysis to describe discrete conformational states of the binding site in which ligands could explore different binding modes. The Glide program was used to perform an exploratory ligand docking and, in order to establish a general approach, we selected the Glide standard precision (SP) protocol.57–59 The described approach was initially applied for the re-docking of the two THS ligands that bind the mAhR, in order to test its ability in reproducing the reference binding modes derived from the HIF-2α templates and relaxed in the mAhR binding site (Figure 3). The starting conformation of each THS ligand is the global minimum conformer (obtained as described in the Experimental Section); in both cases it belongs to the same conformational cluster of the crystallographic structure. To broadly explore the conformational variability of the flexible THS ligands, the first 10 docking poses for each representative mAhR conformation were analyzed. Two clusters of poses were obtained for both ligands and the most populated cluster (60% for THS-017, 73% for the THS-020) reproduces the reference binding mode. The pose with the lowest RMSD to the reference obtained in each representative mAhR conformation was retained and subjected to energy minimization, according to the protocol described in the Experimental Section. Finally, the energy minimized complexes were ranked on the basis of the RMSD to the reference. The best final poses (Figure SI3) reproduce well the ones obtained from the templates, with global RMSD values of 2.74 Å, for the THS-017, and 0.43 Å, for THS-020. Only a conformational rearrangement of the THS-017 thiophene ring, due to the flexibility of the connecting methylene group, is observed. The network of stabilizing interactions with the mAhR previously discussed (Figure 3) was maintained in these poses. TCDD is the highest affinity and most potent AhR ligand.6,40 Considering its extremely high potency and toxicity in many species, it has been widely studied in order to assess the structural determinants that drive its strong association with the AhR. The TCDD bestranked docking poses in the representative mAhR conformations, obtained applying our docking protocol, are shown in Figure 4a. Three clusters of TCDD placements within the cavity are observed: an “internal” pose (in mAhR_HM01), two “central” poses (in mAhR_HM70 and mAhR_HM77) and an “external” pose (in mAhR_HM79). Energy minimization of the complexes (see Experimental Section) permitted the relaxation of the TCDD molecule, its 5 Å shell of surrounding residues, and the residues that we experimentally identified as the TCDD binding fingerprint (Thr 283, His 285, Phe 289, Tyr 316, Ile 319, Phe 345 and Ala 375)17. After the energy minimization, almost all the docking Motto et al. Page 7 J Chem Inf Model. Author manuscript; available in PMC 2012 November 28. N IH PA Athor M anscript N IH PA Athor M anscript N IH PA Athor M anscript poses converge in one conformational cluster that occupies the middle of the cavity, suggesting a favorable intermolecular interaction network in this zone. To better investigate the observed positional convergence of TCDD in the representative mAhR homology models, we carried out combined Monte Carlo/Stochastic Dynamics (MC/ SD) simulations.60 The MC/SD procedure differs from a normal dynamics simulation in that it uses a mixture of Metropolis Monte Carlo61 and dynamics steps in order to greatly increase the rate at which a simulation explores conformational space. The “internal” (mAhR_HM01+TCDD) and “external” (mAhR_HM79+TCDD) poses obtained by molecular docking were selected as starting conformations for 1.5 ns of MC/SD simulations. These simulations revealed that TCDD could move within the cavity by 0.4–4.0 Å of the RMSD, relative to its starting conformations. Almost 30% of the TCDD conformational population sampled during the MC/SD showed a hydrogen bond with the Thr 283 side chain and almost 70% forms a hydrogen bond with the Gln 377 side chain, thus confirming a certain degree of stabilization of the ligand in the middle part of the cavity. MC/SD results appear to correlate with both the observed docking poses and the energy minimization results and indicate that TCDD, a highly symmetric and flat molecule, can explore the cavity extension but seems better stabilized within the middle of the binding site where its oxygen atoms can form hydrogen bonds with Thr 283 and Gln 377 side chains. The energy minimized mAhR/TCDD complexes were ranked using Glide extra precision (XP) scoring function53 in order to select a representative complex. The best scoring complex (mAhR_HM01+TCDD) is shown in Figure 4b. The mAhR binding site was mapped using SiteMap program62 in order to verify matching between the ligand and protein properties. An extended hydrophobic channel was identified at −0.5 kcal/mol together with smaller hydrogen bond donor and acceptor globular volumes at −8 kcal/mol. TCDD fills the hydrophobic channel and orients one oxygen atom near the hydrogen bond acceptor spot (Figure 4c). Thus, the observed ligand binding placement seems to satisfy the mAhR physico-chemical requirements. Superimposition of the best scored complexes for TCDD and the THS ligands allowed direct comparison of the binding mode for all the ligands (Figure 5). TCDD and the THS ligands occupy almost the same region in the binding site. However TCDD is projected slightly more into the cavity than the THS ligands. In fact the former is a rigid and flat molecule while the latter ligands assume a stable bent conformation driven by the intramolecular hydrogen bond within the amino and nitro groups. Computational and experimental data correlation analysis The Glide XP scores for the best mAhR/ligand complexes, both obtained by docking and derived from the templates, are summarized in Table 5. In both cases, the computational ranking is coherent with the differences observed in the experimental ligand binding affinities between the TCDD and the THS ligands. However, the scoring function applied is not sensitive enough to clearly differentiate the higher intermolecular stabilization of THS-017 with respect to THS-020 in the binding site. Extensive experimental mutagenesis of the mAhR LBD has been carried out over the last few years in order to identify and characterize the key residues involved in TCDD binding and these experimental data were used to assess the binding pose of TCDD obtained by molecular docking. We focused our attention specifically on the residues contained within the 5 Å shell around the docked ligand (Table 6) and indeed we found that these include all the residues previously identified as the TCDD binding-fingerprint (Figure 4b), for which mutagenesis resulted in the complete loss of TCDD binding and TCDD-induced DNA binding.17 Other interesting residues that lie within the 5 Å shell around the docked TCDD were identified by site-directed mutagenesis and functional analysis, carried out by our Motto et al. Page 8 J Chem Inf Model. Author manuscript; available in PMC 2012 November 28. N IH PA Athor M anscript N IH PA Athor M anscript N IH PA Athor M anscript group and others, as important for TCDD binding. In particular: mutation of the residues Pro 291, Phe 318 and Ala 361 reduced or eliminated TCDD binding;17,63,64 mutation of Cys 327 to Ala partially reduced TCDD binding;17 Met 334 and Gln 377 mutagenesis had variable effects, with Q377L63 and M334A17 having slightly or partially reduced TCDD binding, Q377A with a significantly greater reduction in binding, and M334E with no ligand binding.17 Mutagenesis of Phe 281 did not appear to affect TCDD binding.63 Analysis of the new mAhR LBD homology model and of the TCDD docking also identified other interesting residues that have not been previously examined by mutagenesis and functional analysis (Table 6). Several of these residues (Leu 302, Leu 309, Leu 347 and Ser 359) were selected for experimental mutagenesis analysis as they project their side chains toward the two sides of the elongated binding cavity, with the three leucine residues lining the inner cavity and Ser 359 present at the entrance of the cavity. Accordingly, we generated a series of mutant AhRs containing individual alanine substitutions at these specific positions (L302A, L309A, L347A and S359A) and demonstrated that these specific mutations did not negatively affect expression levels of the AhRs as each had levels of the in vitro expressed protein similar to that of wtAhR (Figure 6a). All the mutations produced significant effects on TCDD binding: AhRs containing L302A, L309A or L347A mutations exhibited loss of [3H]TCDD binding or dramatically reduced ligand binding activity (to less than 50% of the wtAhR), and ligand binding to the AhR containing the S359A substitution was reduced by ~40% (Figure 6b). Interestingly, overall reductions in the amount of TCDDinducible AhR DNA binding observed with each mutant AhR were very similar to the percent loss of ligand binding observed with each (compare Figure 6b and 6c). These results strongly suggest that the loss of ligand-dependent transformation and DNA binding is due to loss of ligand binding activity and is not due to alterations in other steps in ligand-dependent AhR transformation (i.e., hsp90 interactions, ARNT dimerization and/or DNA binding). Therefore, these new mutagenesis data further confirm the proposed binding mode of TCDD in the mAhR binding site, particularly highlighting the role of the hydrophobic environment conferred by the three leucine residues in the inner part of the cavity in ligand stabilization. Further validation of the reliability of the proposed approach was obtained by docking TCDD into the human AhR (huAhR) model. We previously highlighted14,65 the critical role of the unique internal residue of huAhR different from a residue of the mAhR, V381 instead of A375 in mAhR (Figure SI1), in determining the lower TCDD binding affinity observed for the huAhR.33 We proposed that the bigger valine side-chain reduces the internal space required for accommodating the TCDD molecule, and confirmed this hypothesis by observing the dramatically reduced TCDD binding to the A375V mutant mAhR.14,65 These findings were also confirmed by molecular docking, carried out by our group and others, into the homology model of this mutant17,18 and of the huAhR.18 The huAhR homology model was built using the three holo HIF-2α structures as templates (alignment in Figure SI1) and four representative conformations were selected from the 100 models generated by MODELLER (selection funnel reported in Scheme 1). The same ensemble docking and minimization procedure adopted for predicting the TCDD binding pose in the mAhR was used (see Experimental Section). The comparison between the Glide XP best scoring pose of TCDD in the huAhR and in the mAhR shows that, even though the modeled cavity volumes were similar (803 Å3 in huAhR and 813 Å3 in mAhR), there are differences in the ligand orientation. In particular, the molecular plane of the TCDD is slightly rotated and shifted towards the entrance of the cavity in the huAhR, compared to that in the mAhR, and this shift is due to the sterical hindrance of the Val 381 side-chain of the huAhR (results not shown). As a consequence, in the huAhR there are increased distances from the ligand to the side-chains for which highly stabilizing interactions were predicted from the TCDD/mAhR pose (in particular Gln 383 and Thr 289, corresponding to Motto et al. Page 9 J Chem Inf Model. Author manuscript; available in PMC 2012 November 28. N IH PA Athor M anscript N IH PA Athor M anscript N IH PA Athor M anscript Gln 377 and Thr 283 in mAhR). Accordingly, the best Glide XP score for the huAhR pose (−7.48) was higher (less favourable) than for the mAhR one (−8.70, in Table 5). Together, these results confirmed the ability of the proposed modeling-docking protocol to correctly rationalize the binding affinity differences for the same ligand within the LBD of the AhRs from different species. Virtual ligand screening of PCDDs To verify the improvement resulting from the use of AhR LBD models derived from holo HIF-2α template structures in ligand docking, as well as to validate the proposed docking protocol, the virtual screening of a set of ligands with known experimental binding affinity data was performed. PCDDs were selected as the ligands for these studies, as they are a class of environmental contaminants that are known to bind the AhR with a wide range of affinities, from the highest known affinity ligand 2,3,7,8-TCDD to inactive compounds, depending on the chlorine substitution pattern on the aromatic rings.5–7 In particular, a group of 14 PCDDs for which a homogeneous set of binding affinity data for the rat AhR is available34,35 were selected for these studies. Given that the best available set of ligand binding data was derived using rtAhR, a rtAhR LBD homology model was developed using the same protocol described for the mAhR and the huAhR models (the reference sequence alignment is reported in Figure SI1). Four representative models were selected and used for the ensemble docking of PCDDs and the following refinement and rescoring steps (see Experimental Section). Moreover, the same procedure was also used to analyze PCDD binding to the rtAhR model we developed starting from the apo NMR structures of the HIF-2α and ARNT PAS B,15,16 for comparative purposes. Similarly to the differences observed in the cavity space of the mAhR apo and holo models (Table 3), the obtained rtAhR binding cavities also have very different CASTp volumes, ranging from 344 Å3 to 773 Å3 for the four representative apo models and from 636 Å3 to 1071 Å3 for the four representative holo models. The Glide XP best scoring poses obtained at the end of the docking, minimization and rescoring procedure for the 14 PCDDs into the apo and holo models of the rtAhR are shown in Figure SI4. It can be observed that in the apo model (Figure SI4 a) favorable poses into the binding cavity were only found for a limited number of ligands, whereas others were excluded, irrespective of their experimentally established binding affinity. In particular, the poses of the lowest chlorinated PCDDs (1-MCDD and 2,8-DCDD) were predicted in the center of the cavity, while those of six other PCDDs (1,2,3,7,8-PeCDD, 1,3,7,8-TCDD, 2,3,6-TrCDD, 2,3,6,7-TCDD, 2,3,7TrCDD and 2,3,7,8-TCDD) were placed at one side of the cavity. The hierarchical series of filters used by the Glide program57,58,66 to search for possible locations of the ligand excluded the other PCDDs from the binding site, probably due to their sterical hindrance. This group includes OCDD and all the PCDDs that have a bulky substituted aromatic ring due to chlorination at both the ortho positions relative to the oxygen atoms (1,2,3,4-TCDD, 1,2,3,4,7-PeCDD, 1,2,3,4,7,8-HxCDD, 1,2,4-TrCDD, 1,2,4,7,8-PeCDD). In contrast, binding poses were found in the cavity of the holo model for all the PCDDs, except OCDD (Figure SI4 b). As expected from the observation that the rtAhR and mAhR modeled binding cavities have similar volumes and share the same internal residues,17 the predicted 2,3,7,8-TCDD pose lies in the same central region as in the mAhR model and interacts with the same surrounding residues, whose role in TCDD binding to the mAhR was validated by mutagenesis analysis. The other PCDDs were docked in the same zone but with different orientations of the molecular planes and intermolecular interactions; in particular, the poses of the PCDDs that were excluded from the cavity of the apo model were placed with the bulky aromatic ring close to the large entrance of the holo cavity. It is conceivable Motto et al. Page 10 J Chem Inf Model. Author manuscript; available in PMC 2012 November 28. N IH PA Athor M anscript N IH PA Athor M anscript N IH PA Athor M anscript that the successful results reported in docking of ligands into the reduced cavity of AhR models derived from the HIF-2α apo structure18,21 is due to the use of different docking protocols having more tolerant intermolecular potentials or explicit treatment of the receptor flexibility. It is known that docking scoring functions suffer limitations in predicting the binding energies,59 and a rescoring procedure is often required to obtain a better ranking of the binding poses. Indeed, the Glide XP scoring function was not able to correctly evaluate the observed differences in the PCDD poses in the holo rtAhR model and failed to reproduce the experimentally determined ranking of their binding affinities. Therefore, a rescoring of the obtained poses was performed with the molecular mechanics generalized Born/surface area (MM-GBSA) protocol67 (see Experimental Section). Such calculations are largely used in lead optimization for rescoring poses unfavorably ranked by docking programs,21,68,69 and this method has been proven to perform best for congeneric series of ligands.21 The resulting correlation between the MM-GBSA binding free energy and the experimental pEC50 values34,35 for the PCDD set (Figure 7) was very good, with a correlation coefficient R2 = 0.81. Lower correlations with the experimental binding affinities were obtained by other authors, by using similar docking and rescoring programs, for the docking of PCDDs and other classical AhR ligands to an AhR LBD model derived from the apo HIF-2α template.21